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We examine a system of doubly degenerate p-orbital polarized fermions on a two-dimensional 
square lattice with a strong on-site interaction. We consider the system density at the half filling 
limit and tackle the strong attractive interaction using a perturbation theory. We treat the four-site 
square plaquette interaction term generated from the directional tunneling dependence of p-orbitals 
using the fourth order in perturbation theory. We map the strong coupling particle Hamiltonian into 
an effective spin-Hamiltonian and then use a variational mean field approach and a linear spin-wave 
theory to study the phase diagram. Further, we discuss the experimental signatures of these phases 
within the context of current cold-atom experimental techniques. 


I. INTRODUCTION 


Cold atoms loaded into an optical lattice formed by interference of counter propagating laser beams are regarded 
as a quantum simulator for many-body condensed matter systems [lj. The flexibility of generating various lattice 
geometries and the controllability of the dimension, lattice parameters, and interaction parameters within the cold 
atomic setups make them ideal test beds for studying strongly correlated many-body phenomena. In condensed matter 
many-body systems, such as transition metal oxides and rare-earth materials, both spin and charge degrees of freedom 
play an important role determining their properties. For some materials, the orbital degrees of freedom is also an 
vital factor. The orbital degrees of freedom is active when the electronic orbitals are degenerate and partially filled 
with electrons. The strongly correlated nature and spatial anisotropy of d and / orbital electrons in some correlated 
materials are responsible for various exotic magnetic, superconducting, and transport properties. 

Due to the recent advancement of laser technology and the progress in experimental techniques, cold atomic physics 
and condensed matter physics are in close proximity to each other. Condensed matter models can be engineered using 
atoms trapped in a combined harmonic and lattice potentials using various laser arrangements and other experimental 
laser techniques in cold atom setups. These so-called quantum simulators or atoms in optical lattices, not only provide 
better opportunities to explore existing condensed matter many-body phenomena, but also allow one to seek and 
explore new states of matter and their physical properties. Unlike solid state electronic models, the optical lattice 
setups are rigid and robust. As a result, the Jahn-Teller distortion and the conventional phonons are absent. However, 
these effects can be engineered using various laser techniques and different types of atoms at low temperatures. The 
electrons which are fermions are the constitute particles in condensed matter lattice models. However, not only 
fermions but also bosons, as well as mixture of fermions and bosons can be used as constitute particles in optical 
lattice models. Opposed to the condensed matter systems, even p-orbital occupation of atoms in optical lattices 
can show strongly correlated many-body phenomena. Recently, a tremendous experimental and theoretical progress 
has been achieved in the direction of understanding p-orbital bands of both fermions and bosons in optical lattices. 
Fermions on a p-orbital band have been imaged and probed by transferring them into the p-orbital band using a 
sweep across the Feshbacli resonance Q. As we consider in this letter, p-orbital bands can be activated even simpler 
way by completely filling the lowest s-orbital band with two pseudo-spin 1/2 atoms and then having more than 
two atoms per site. The population of p-orbital band bosons were observed in moving lattice experiments [3] and 
dynamically deformed double-well lattices Q. The p-orbital boson decay time has also been measured by pumping 
atoms by Raman transitions Q . The quantum mechanical phenomenon associates with orbital degrees of freedom of 
atoms in optical lattices can be captured by a p-orbital band Hubbard Hamiltonian j6b& . The emergence of various 
quantum phases of p-orbital band bosons in optical lattices has been recently proposed |9H27j] . A tremendous progress 
in studying orbital exchange driven magnetic properties of p-orbital band fermions has also been reported [28l - l43l |. 

In this letter, we study the phase diagram of p-orbital attractive polarized fermions loaded into a optical lattice. 
The p-orbital occupation can be achieved by having two one-component fermions at each site. In this scenario, the s- 
orbital is completely filled and the s-orbital atom remains inert at the low energy subspace. The attractive interaction 
can be generated simply by adjusting the two-body scattering length to be negative values using Feslibach resonance. 
We consider the system density at the half filling and strongly interacting limits. At the half filling density limit, on 
average, only one atom resides in a p-orbital at each site. First, we map the system Hamiltonian into an effective spin 
Hamiltonian. The fourth order tunneling process favored by the spatial anisotropy of the p-orbital orientation is taken 
into account using the fourth order in our perturbation theory. Second, we map out the magnetic phase diagram of 
the effective spin Hamiltonian using a variational mean field theory. Third, we use a linear spin-wave theory to study 
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the phase diagram and magnon excitation spectrum. 

The letter is organized as follows. In section II, we introduce the effective model and discuss the details of mapping 
it into an effective spin-model. In section III, we discuss our variational mean field approach and provide the resulting 
phase diagram of the effective spin-model. In section IV, we provide the details of our linear spin-wave theory for the 
effective spin-model and discuss the phase diagram and magnon excitation spectrum. In section V, we discuss the 
connection of the resulting phase diagram to the original model and experimental signatures of the magnetic phases. 
Finally in section V, we draw our conclusions. 


II. EFFECTIVE MODEL 

We consider multi-orbital system of polarized fermions on a two-dimensional square lattice. We assume that the 
atoms are loaded to the lattice such that the s-orbital is completely filled and the dynamic of the system is determined 
only by the p-orbital atoms. Further, we assume that only p x and p y orbitals are active on the two-dimensional (2D) 
lattice. This scenario is reasonable as the s and p-orbital bands are separated by a large band gap at experimentally 
relevant parameters and p z orbital has a higher energy in 2D configurations. By imposing a strong laser beam along 
the ^-direction, the p 2 -orbital band energy can further be increased. As sketched in FIG. [U the p x orbitals overlap 
only in the x-direction. As a result, the tunneling between p^-orbitals in ^-direction is much greater than that of the 
y-direction. Similarly, the hopping between p y orbitals along y-direction is greater than that of the ^-direction. By 
representing the p x orbital occupation by pseudo-spin f and p y orbital occupation by pseudo-spin 4, the dynamics of 
the spin polarized atoms in pseudo-spin-dependent 2D optical lattice can be described by an effective Hamiltonian, 


H = - “ U Yl( ni t ~ V 2 )( n 4 - V 2 ) ~ 4r C i<r + V Y2 n i n j- (1) 

<ij>,cr i icy <.ij> 

The first term is the kinetic energy and is proportional to the tunneling amplitude t l J along the bond ij = x,y 
directions of a =f {p x ),i ( p y ) atoms. The operator c\ a (ci a ) creates(destroys) a Fermi atom with pseudo-spin a =t,4- 
at lattice site i. The second term describes the on-site interaction energy U > 0. The density operator or the 
occupation number operator is rn a = c\ a Ci a and p is the chemical potential. Notice < ij > indicates only the nearest 
neighbor pair of sites and we neglect tunneling beyond the nearest neighbors. The last term represents the nearest- 
neighbor weak, off-site repulsive interaction. All other nearest neighbor interactions, such as bond-charge, spin-spin, 
and pair-hopping are neglected. 

We consider the system at the half filling density, so that on average, only one atom resides at each lattice site 
occupying either p x or p y orbital. Due to the directional dependence of p-orbitals, we take the tunneling amplitude 
of atoms t l J as t l J~' x = t l ^~ v = t and tV~ x = t l J~ v = at with 1. As shown by the Hamiltonian above, we assume 
on-site interaction is attractive (—U < 0), which can be achieved by tuning the free-space two-body scattering length 
to be negative values. The directional pseudo-spin tunneling amplitudes and the on-site attractive interaction favor 
a four-site ring-exchange or plaquette interaction term in the form, 


Hp — K c ia C Jcr ,c l a 'Qat 

<ijkl>,cr 


( 2 ) 


where < ijkl > represents the four-site square plaquette and the four spin coupling constant K oc t 4 /U 3 can be 
derived in fourth order hopping process. This process is schematically shown in FIG. [2j 


Psudo-spin Model at the strong coupling Limit 

In the strongly interacting limit where t <C V < U , the tunneling term acts as a small perturbation. In this section, 
we derive an effective pseudo-spin Hamiltonian in fourth order perturbation theory. The local Hilbert space has four 
states, namely empty site denoted by |e), singly occupied sites with up | f) or down | 4-) spins, and a doubly occupied 
site with opposite pseudo-spins \d). The energies of these states in unperturbed Hamiltonian, i.e the on-site part of 
the Hamiltonian, are given by Ei = ~U/ 4, Ei = —p + U/ 4, Ei = —p + 17/4, and Ej = —2 p —U/ 4, respectively. Note 
that our notations up and down spins represent the p x and p y orbital occupation. At the half filling density limit 
where p = 0, the singly occupied sites have higher energies and the empty and double occupied sites energies are 


3 



FIG. 1: While the tunneling amplitude t along the x-direction between the p x orbitals is large, the tunneling 
amplitude at, (with a -C 1) between the p y orbitals is much weaker. The gray region represents the p x orbital, 

whereas the white region represents the p y orbital. 
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FIG. 2: Four site plaquette interaction at the half filling density limit for strongly attractive fermions due to the 
spatial anisotropy of p orbitals. While the up-arrow (t) represents the p x orbital occupation, the down-arrow (4_) 

represents p y orbital occupation. 


degenerate. The local Hilbert space of two-neighboring sites has ten states, however six of those states that involve 
a single site have higher energies. Therefore, the low-energy subspace constitutes only the four states, |e, e), |e,d), 

| d, e), and | d, d) with energies E = —U/ 2, E = —2/z — U/2, E = —2/z — U /2, and E = — 4/z — U/2 + 4V, respectively. 
Here | e,d) stands for two neighboring sites with a empty and a doubly occupied sites and other notations have the 
similar meaning. 

In the second order perturbation theory, the energy correction due to the tunneling between two neighboring sites 
can be represented as a matrix M, 


/ ( e > e l \ 

1 ( 

{e,d\ 


(d,e j 


\(d,d\) 

V 


-U/2 

0 

0 

0 


0 

-2/z - U/2 - (a 2 + 1 )t 2 /{U + V ) 
—2at 2 /(U + V) 

0 


0 

-2 at 2 /{U + V) 

—2/z - U/2 - (a 2 + l)t 2 /(U + V) 
0 


0 \ 

/ l e i e ) \ 

0 

1 e,d) 

0 

\d, e ) 

4/z - U/2 / 

\\d,d)J 


Here we neglect the fourth order corrections as the dominant fourth order correction is taken into account through 
Hamiltonian Hp. The fourth order correction here is much smaller as it is proportional to a 2 or its higher powers. 
The T in the column vector in Eq. (3) indicates the transpose. At the half filling density limit, we map this low- 
energy Hamiltonian into an effective spin-1/2 system by mapping the local double occupied sites by pseodo-spin-ups 
and the local empty sites by pseodo-spin-downs, ie |d) —> | f)s P m and |e) —> | 4-}spmj such that the z-component 
of the spin operators S z = ha z /2 acts as a z \ t) sp m = +| t)s P m and a z \ 4.} sp i n = —| i) sp in■ Notice that these 
new pseodo-spin states | t)s P m and | 4-}s P m are different from our previous pseodo-spin states | f) or I 4-} that 
represented the p-orbital occupations. We start with the most general two site effective Hamiltonian Hij that respects 
the global SU(2) symmetry. Up to the quadratic terms, the spin-1/2 effective Hamiltonian must have the form 
Hij = J2 U {A V S/Sj + B^S/} + C, where v = x, y, z. By calculating the matrix elements with respect to this effective 
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Hamiltonian Hij and equating them to the matrix elements in Eq. (3), all exchange interactions terms, A v , B u , and C 
are derived. In terms of Pauli matrices a that represent the new pseodo-spins, the effective pseodo-spin Hamiltonian 
is then written as, 


H = 


<ij> 


{-J(afaJ + afaJ) + A<a~} + | £ a*a*a* k af, 

<ijkl> 


(4) 


where J = at 2 /{U + V), A = (a 2 + l)t 2 /[2(U + V)] and K = 12t 4 U/[(12V + 3U) 3 (8V + 3C/)]. The last term is 
originated from the fourth order tunneling process shown in FIG. [2] All other fourth order tunneling processes are 
much weaker for the present system. We have already set the half filling condition p, = 0 and ignored the constant 
energy shift of C = —U/2 — (a 2 + l)t 2 /[2(U + V)]. At half filling all the linear terms B v vanish. 

Without the last quartic term, Eq. (4) represents the well-known ferromagnetic XXZ model. As the model contains 
both ferromagnetic exchange J and antiferromagnetic exchange A, the spin Hamiltonian is frustrated even on a 
geometrically non-frustrated two-dimensional square lattice with only nearest neighbor exchange interactions. This 
model can show finite temperature Berezinskii-Kosterlitz-Thouless transition without finite magnetic order and a 
Ising type phase with long range order 44]. The XXZ spin model is equivalent to the hard-core boson Hamiltonian 
with nearest neighbor hopping t = 2J and nearest neighbor repulsive interaction V = A [ihj . In bosonic language, 
the magnetic order in the XY plane corresponds to superfluid order, while the magnetic order in the Z direction 
corresponds to density order. 


III. A VARIATIONAL MEAN FIELD THEORY FOR THE SPIN HAMILTONIAN 

In this section, we use a variational mean field approach to study the phase diagram originating from the competition 
between three terms in the spin Hamiltonian given in Eq. (4). Notice all coupling constants, J, A, and K are positive. 
While J term favors a ferromagnetic spin ordering inj-j/ plane in spin space, A term prefers an antiferromagnetic 
ordering in z-direction. However, the K term discourages both of these ordering. 

First, we break up the original square lattice into two interpenetrating square sublattices A and B. Next, we take 
our normalized variational density matrices for sublattice A and B as 

Pa,b = ^ + ^m(smda x ± cos 9a z ), (5) 

where m and 9 are variational parameters and the upper sign is for the sublattice A and the lower sign is for the 
sublattice B. This choice gives us the sublattice magnetization itia,b = Tr(pA,B&) = Fmcos9z + msin0i, here 
Tr(L) represents a trace of a 2 x 2 matrix L. Thus, m = 0 represents the paramagnetic state and m / 0 represents 
three magnetically ordered states depending on the value of 9. For 9 = 0, the system is an antiferromagnetic state 
with ifiA = —friB and its staggered magnetization is in the z-direction. For 9 = 7 t /2, the system is a ferromagnetic 
state with itia = ms and its magnetization lying along the azdirection. On the other hand when 9 ^ 0 , the system 
is a canted antiferromagnet with uia ^ in b- 

The free energy F = E — ksTS with energy E and entropy S, can easily be calculated using E = Tr(pH) and 
S = —Tr(p In p), where ks is the Boltzman constant and T is the temperature. Defining dimensionless free energy per 
site f(m, 9) = F/[JN\ with energy scale J = t 2 /[U+V] and N being the number of lattice sites in the two-dimensional 
square lattice, we find 


/(to, 9) = —2 am 2 



1 + TO 
2 


In 


1 + TO^ 

) 



( 6 ) 


Here the dimensionless temperature r = ksT/J and the dimensionless four-site square plaquette coupling n = K/ J. 
The minimization of the free energy with respect to two variational parameters, df /dm = 0 and df /d9 = 0 leads to 
two coupled mean field equations, 


— 4am — 2(a 


1) 2 to cos 2 9 + Km 3 cos 4 9 + 



1 + TO 


= 0 


1 — TO 


(7) 
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(a) The phase diagram at a fixed dimensionless temperature r = 1. 


(b) The phase diagram at a fixed tunneling anisotropy a = 0.3 


FIG. 3: Magnetic phase diagram obtained from the variational mean field theory of Hamiltonian (4). Theory 
predicts Antiferromagnetic (AFM), canted antiferromagnetic (C-AFM), and paramagnetic (PM) phases in a — n 
t — k planes, where a is the tunneling anisotropy and k is the the dimensionless four-site square exchange 

interaction parameter. 


and 


— to 2 sin 6 cos 6[nm 2 cos 2 9 — 2(a — l) 2 ] = 0. (8) 

By simultaneously solving these mean field equations, we find 


f 0, if (a — l) 2 < 0; 

cos26,= l 2( °J )2 > if 0 < (a — l) 2 < nm 2 /2; (9) 

[ 1, if (a — l) 2 > Km 2 / 2. 

Obviously, the first condition cannot be fulfilled, thus a a;y-ferromagnetic phase is not stable for our effective spin- 
Hamiltonian in Eq. (4). Using the theta values in Eq. (0) and the first mean field equation, we find a canted 
antiferromagnetic (C-AFM) state for k > 2(1 — a) 2 and r < t c where the critical dimensionless temperature is given 
by, 


= wa q-oiv* _ (10) 

tank 1 a / 2(1 — a) 2 / k 

On the other hand, for k < 2(1 — a) 2 , we find an antiferromagnetic (AFM) state if r < 2(a 2 + 1). If any of these 
conditions are not satisfied, the system remains in a paramagnetic (PM) state. The phase diagrams resulting from 
these conditions for various parameter regions are shown in FIG. [3] 

IV. A LINEAR SPIN-WAVE THEORY FOR THE SPIN HAMILTONIAN 

In this section, we investigate the quantum phase transition between AFM state and C-AFM state using a linear 
spin-wave theory (LSWT). In quantum spin systems, the low energy or long-wave length magnetic excitations spectrum 
is formed by spin-waves or magnons. Magnons are collective modes associated with the coherent precession of spins 
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at a given site. The original introduction of spin-wave concept was proposed by Bloch and Slater 46]. Holstein and 
Primakoff later introduced the quantization of spin waves in terms of bosonic operators that represent magnons |47 |. 
LSWT is the linearized version of this bosonic theory and it is a well known tool in the study of quantum spin 
systems. In this method, each spin operator is represented by two Holstein-Primakoff (HP) bosons 48}. The vacuum 
state of Holstein-Primakoff bosons will be a broken symmetric state, therefore the theory is valid only for long-range 
magnetically ordered states. Some generalizations have been proposed to validate LSWT for short-range magnetically 
disordered states 49|-52j]. Here we assume that the spin system is in either AFM ordered state or C-AFM ordered 
state and study the phase transition between them. In LSWT, the quantum spin Hamiltonian is mapped onto an 
interacting boson Hamiltonian and then keep only the quadratic terms by neglecting the higher order bosonic terms. 
The higher order bosonic terms represent the interaction between HP bosons or magnons. Consequently, LSWT is a 
method of converting an interacting spin Hamiltonian into a non-interacting bosonic Hamiltonian. 

First, we rotate the spin operators a® into the local reference frame at each site using Euler rotation cq = U_i&i 
[here <r°’s are the original spin operators in Eq. (4) and oVs are the spin operators in local reference frame], where 
the rotation matrix V_ i: is given by 53, 


Hi = 


0 1 0 

sin 6i 0 cos 9t 


( 11 ) 


Notice that we consider only the in-plane rotation so that we have set the second Euler angle xp to be zero. In the 
local frame, the Hamiltonian in Eq. (4) reads H = jeB) Hab + H q , 


HaB — jx.cAx&Bx + Jy& Ay&By + Jz&Az&Bz + Jzx&Az& Bx + Jxz&Ax&Bz, 


( 12 ) 


where the four spin interaction terms H q = J q &® z v 0 j z a kz a iz with J q = A'/4, in the local frame is given by, 


Hq = J q ^2 l cos 2 9A cos 2 9 B &i Z <jj Z <Jkz&iz - 2 sin(2 6 A ) cos 2 9 B [cri Z <Jj Z a kx ai z + cr ix aj z a k zcriz\ (13) 

ieA,jeB,keA,ieB 1 

+ sin 2 9 A cos 2 9 B a ix cr jz a kx ai z - ^ cos 2 9 A sin(2 9 B )[a iz a jz a kz ai x + d iz a jx a kz ai z } 

+ i sin(26bi) sin(20 s ) [ a iz &jz&kx&lx + &ix&jz&kz&lx H - &iz&jx^kx^lz “1“ &ix&jx&kz&lz\ 

“ sin Qa sin(2$£) \o‘i x CTjz^kx^’lx ^ix^jx^kx^lz] COS Oa sin OB^iz^jx^kz^kz^lx 
-- sin(2 6 A ) sin 2 0 B \(*iz Gjx&kx@lx “1“ &ix &jx&kz&ix] T sin - 9a sin” 9 B (Ti X (jj X (j kx <Ji x 1 . 


The anisotropic coupling constants J x = — J cos 9a cos 9 b + A sin 9a sin 9 B , J y = — J, J z = — J sin 9a sin 9 B + 
Acos9 a cos 9 b , J zx = — J sin 6 a cosd B — A cos 9a sin 9 B , and J xz = — J cos 9 a sin 9 B — A sin 9a cos 9 b . Notice that we 
have divided the two-dimensional lattice into even and odd sublattices, denoted by A and B such that 0* = 9a and 
9j = 9 b for i € A and j £ B, respectively. 

Second, we introduce Holstein-Primakoff bosonic operators a* and bj for sublattices i € A and j € B via Holstein- 
Primakoff transformation, 


for the sublattice A and, 


Viz = cr ~ a] cii and cr + 



Tj z = bfij - a and ct+ 



(14) 


(15) 


for the sublattice B, where n* = aja,; and rij = b^bj are the boson number operators and eq + = a x -(- ia y = (a l )1 is 
the spin raising operator. 
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Third, we express the spin Hamiltonian in terms of HP bosons and keep only the terms up to the quadratic order 
in bosonic operators. Finally, by introducing the Fourier transformation of the bosonic operators, 



where N is the total number of lattice sites, we write the linear spin-wave Hamiltonian as H = Hq + + H 2 . The 

zeroth order that describes the classical energy is given by, 


s 2 Nz s A Nz 

Hq = — [— J sin 9a sin 9 b + A cos 9a cos 9b] d- —— J q cos 2 9a cos 2 i 


(17) 


where z = 4 is the number of nearest neighbors and s = 1/2 is the effective spin. By minimizing the classical energy 
with respect to 9a and 9b, we find two relations, 


2 J xz + J q S 2 sin(26U) cos 2 9 B = 0 (18) 

2 J zx + J q S 2 sin(20s) cos 2 9a = 0. 

By solving these equations simultaneously, we find two solutions 9a = 9b = 0 or 9a=9b = cos repre¬ 
senting the AFM state and C-AFM state respectively. This indicates that the spin system undergoes a second order 
phase transition from an AFM state to a C-AFM state as one increases J q . The phase boundary between these two 
states is given by the condition J + A = J q s 2 and the spin canted angle with respect to AFM spin orientation is given 

by 9 C = cos-M^]. As a demonstration, we plot the classical energy per site as a function of canted angle for two 
different sets of representative parameters in FIG. [I] As can be seen, the classical energy shows a local minimum at 
zero canted angle for smaller values of J q , where it shows a local minimum at a finite value of canted angle ( 9 C ) for 
larger values of J q . 

The first order terms that are linear in creation and annihilation bosonic operators are given by, 

E E E E («-*+4) as) 

ij k ij k 

J S ' 7 / 2 

+ 1 ^i /2 cos2 sin(20s) ££^' (b k + bl k ) 

ijlm k 

— cos2 sin(26(4) ££e-^(a- fc + at). 

ijlm k 

As the spin system is in either AFM or C-AFM ground states with the angles that minimize the classical energy given 
by Eq. (18), all the first order terms given in H\ vanish. 

The second order terms of the spin wave Hamiltonian is given by, 


H 2 = 


£ 


Jo{a\ a k + b{b k ) + Ji(a fe a_ 


k + 4 a -fc 


b k b-k 


^I^-fc) J~( a kb*_ k + a \(b— fc) + J+(dkbk 


A b l) 


( 20 ) 


where the momentum dependent coefficients in each term are, Jo = sJ z z + zs 3 J q cos 2 9b (2 cos 2 9a + 7 fcsin 2 6 bt), 
Ji = zs 3 Jq'yk sin 2 9a cos 2 9b, = [J x s — J y s — s 3 J q sin(26U) sin( 20 s)] 2 ( 7 fc/ 2 , and J + = [. J x s -f J y s — 

s 3 J q sin(20^) sin( 20 s)] 2 ; 7 fe/ 2 . Here the structure factor 7 = 1 /z^2$e lk ' S is introduced with the nearest neighbor 
vector S = dx + dy , where d is the lattice constant. Finally, using the equation of motion for bosonic operators 
e k = {a k ,b- k ,a' k ,b'_ k }, 















Oc 

FIG. 4: Classical energy derived from linear spin wave theory at A = 2 J. While the black solid line is for J q = J 
and the gray solid line is for J q = 25 J. Notice the minimum energy at canted angle 8 C = 0.805 radians for J q = 25 J, 
where as the energy is minimum at 9 C = 0 in the AFM state at J q = J. 


^ — [-^ 2 ? — ^ek&ki ( 21 ) 

we find two distinct eigenvalues A a k and A bk, where [X, F] stands the commutator for operators X and Y. We then 
determine the linear spin wave frequency uik = V^ak^bk, 


u k = //q 2 - (-4 - J-)- (22) 

This magnon dispersion for the AFM state and the C-AFM state is plotted in FIG. [5j We use two sets of representative 
parameters to show the key features in the dispersion relevant to each of these phases. Notice that for some parameters, 
the magnon dispersion can be complex. This complex spin wave frequency indicates the instability of the ordered 
phase and represents the paramagnetic phase. This is not surprising as our spin wave theory assumes the long range 
magnetic order. 


V. CONNECTIONS TO EXPERIMENTS 

Here we considered a two-dimensional optical lattice with fully occupied lower energy s-orbital atoms and partially 
occupied higher energy p-orbital atoms. Such a system can be realized by loading single hyperfine state fermionic 
atoms, such as 40 K or 6 Li. After the s-orbital is fully occupied, partial p-orbital occupation is attainable at a higher 
density filling factors. As the anisotropic nature of the p-orbitals, the Bloch band formed by the p z orbital has a higher 
energy in two dimensions, thus the p x and p y orbitals are degenerate. The strong attractive interactions between 
spinless fermions can be achieved by p-wave Feshbach resonance. 
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FIG. 5: The magnon dispersion Uk for AFM and C-AFM states. 


The paramagnetic (PM), antiferromagnetic (AFM), and the canted-AFM phases elevated by the mapped spin 
Hamiltonian given in Eq. (4) correspond to different particle occupation patterns relevant to the original particle 
Hamiltonian given in Eq. (1). The PM phase represents the random single particle occupation of p x or p v orbitals at 
each site. The AFM phase represents the alternative empty and double occupation, (0, 2, 0, 2,...) pattern. Meantime, 
the C-AFM phase represents the random empty and double particle occupation, (0,2,2, 0,0,..) pattern. The canted 
angle measures the randomness and it depends on the microscopic parameters in the Hamiltonian. In both AFM and 
C-AFM states, every lattice sites are either empty or doubly occupied. Therefore, the predicted magnetic phases of 
the effective spin Hamiltonian can be directly detected by the density pattern in the lattice. These density patterns 
can be directly probed by in situ imaging techniques [55|. Alternatively, one can use single-site detecting techniques 
to map these density patterns experimentally [56|,l57|. At low temperatures, the two Fermi atoms in different orbitals 
at each double occupation sites can bind together to form a boson. As one lower the temperature of the optical 
lattice, these bosons can condensate and show normal-superfluid phase transition. The details of this phase transition 
is beyond the scope of this paper and we leave further details to a separate study. 


VI. CONCLUSIONS 

We have studied the phase diagram of doubly degenerated p-orbital attractive fermions on a square lattice. We 
considered the strongly interacting limit and the half Filing density limit, and derived an effective spin Hamiltonian 
using a perturbation theory. We find that the fourth order square plaquette interaction term generated from the 
anisotropic tunneling process is significant at these limits. We used a variational mean field approach and a linear 
spin-wave approximation to map out the phase diagram of the effective spin Hamiltonian and find three distinct 
magnetic phases. Further, we discussed the experimental connections to these phases and argued that the resulting 
phases can be detected by using currently available experimental techniques in cold gas experiments. 
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